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Summary 

This  report  describes  how  to  invert  seismic  data  collected  with  sonobuoys  in  littoral  inarine 
environments  to  obtain  sea  floor  geo-acoustic  properties.  The  inversion  method  is  an 
adaptation  of  simulated  annealing  (SA)  which  can  rapidly  and  automatically  find  the  best 
earth  model  for  each  data  set.  The  established  SA  algorithm  was  enhanced  with  envelope 
fitting  and  layer  stripping.  The  first  example  is  a  synthetic  data  set  to  demonstrate  (he 
method  for  a  known  earth  model.  The  second  example  is  at  a  750  m  water  depth  site  using 
a  watergun  source  which  has  a  very  clean  and  sh^  pulse.  The  offsets  were  accurately 
determined  prior  to  the  inversion  ciculation  for  this  case  and  a  clear  subbottom  reflector 
was  modeled  in  only  5323  iterations.  The  third  example  was  at  a  site  with  250  m  water 
depth  using  an  airgun  source.  The  direct  airgun  signal  saturated  the  sonobuoys  at  small 
offsets  causing  distortions  in  the  waveforms  and  amplitudes.  Data  offsets  could  not  be 
accurately  measured  and  had  to  be  included  in  the  inversion  problem  which  tripled  the 
number  of  variables.  This  more  difficult  third  example  was  solved  with  35421  iterations. 

This  demonstration  is  with  highly  nonoptimal  data.  The  receivers  are  sparsely  and 
irregularly  spaced  and  have  too  limited  a  dynamic  range  for  the  source  strengths  and  sea 
floor  reflections  encountered.  However,  by  adapting  the  inversion  technique  to  the 
challenges  presented  by  the  field  data,  the  geoacoustic  parameters  of  the  sea  floor  were 
measured  using  rapidly  deployed  inexpensive  and  expendable  equipment.  The  geoacoustic 
properties  that  I  modeled  include  the  sea  floor  depth,  the  water  sound  velocity  and  the 
compressional  wave  velocity  of  the  uppermost  150  m  of  sediment.  In  the  last  case  I  also 
inverted  for  the  sonobuoy  to  airgun  offsets. 

The  reliability  of  the  sonobuoy  -  airgun  and  sonobuoy  -  watergun  data  collection  and 
the  SA  inversion  method  is  demonstrated  with  favorable  comparisons  with  high  quality 
multichannel  seismic  (MCS)  data  from  nearby  sites.  All  of  the  major  Vp  -  depth  features 
that  can  be  expected  to  be  modeled  from  the  field  data  correspond  we|l  with  measurements 
from  the  MCS  data. 


Introduction 

The  purpose  of  the  inversion  technique  presented  in  this  report  is  to  rapidly  obtain  sea  floor 
geoacoustic  parameters  to  support  Navy  ASW  operations.  This  technique  uses  data  that  can 
be  rapidly  collected  from  ships  and  (potentially)  aircraft  using  existing  inexpensive  and 
easily  deployed  equipment.  This  techidque  has  been  previously  demonstrated  using  other 
types  of  marine  seismic  data.  For  example,  in  Wood  and  Lindwall  (1996),  this  technique  is 


applied  to  high  frequency  vertical  acoustic  profiles  of  the  sea  floor,  Lindwall  et  al.  (1995) 
applies  this  technique  to  deep  tow  multichannel  seismic  data  and  Lindwall  (1995) 
demonstrates  this  same  inversion  method  on  the  vertical  acoustic  profiles,  the  deep  tow 
multichannel  data  as  well  as  simulated  seismic  refraction  data. 

The  method  described  in  this  report  is  available  as  a  set  of  FORTRAN  programs  that 
can  analyze  field  data  in  the  proper  (relatively  simple)  format.  The  field  data  must  first 
corrected  for  shot  time  and  shot  to  receiver  offsets  as  well  as  possible.  A  standard  seismic 
processing  software  package  such  as  Seismic  Unix  (available  from  the  Colorado  School  of 
Mines  at  http;//www.cwp.mines.edu/cwpcodes)  is  ideal  for  this  task.  Accurate 
measurements  of  the  water  sound  velocity,  water  depth  and  source  and  receiver  depths 
need  to  be  made  in  the  field  for  a  fast  and  accurate  inversion  solution.  Unknown  or 
inaccurate  environmental  variables  can  be  inverted  for  but  slow  down  the  process.  This  is 
shown  in  my  third  demonstration  where  the  shot  offset  variables  are  included  in  the 
inversion  and  turn  a  30  minute  calculation  into  a  10  hour  calculation. 

The  inversion  is  first  applied  to  synthetic  data  as  a  test  and  then  to  two  sets  of  field 
data.  The  field  data  was  collected  in  the  STRATAFORM  area  near  the  mouth  of  the  Eel 
River  in  California.  The  STRATA  FORmation  on  Margins  program  is  a  coordinated  multi¬ 
investigator  study  of  continental-margin  stratigraphy  initiated  by  the  Office  of  Naval 
Research. 


Methods 

The  rapid  inversion  method  described  here  is  adapted  from  Simulated  Annealing  (SA).  SA 
inversion  uses  synthetic  data  from  a  forward  modeling  algorithm  and  an  evolutionary  test 
criterion  to  determine  the  best  fitting  physical  model  (environmental  parameters)  for  a  given 
data  set.  SA  inversion  is  used  in  cases  where  the  physical  model  cannot  be  directly 
determined  from  the  data  and  the  model  space  is  too  complex  for  the  best  solution  to  be 
found  by  iteratively  moving  from  a  random  starting  solution  to  a  better  one.  The  forward 
algorithm  used  in  this  demonstration  is  a  ray  trace  code  that  uses  a  one  dimensional 
environmental  model,  includes  shear  wave  conversion  losses,  but  does  not  include 
compressional  wave  multiples  or  the  sea  surface  reflections  (Chapman,  1976;  Cerveny  et 
al.,  1977).  It  was  chosen  for  this  demonstration  because  it  is  extremely  fast  relative  to  full 
wavefield  methods  such  as  reflectivity  (Fuchs  and  Muller,  1971). 

The  original  SA  algorithm  was  set  forth  by  Metropolis  et  al.  (1953)  and  reintroduced 
by  Kirkpatrick  et  al.  (1983).  The  first  application  of  SA  inversion  to  velocity  estimations 
from  seismic  waveform  data  is  by  Basu  and  Frazer  (1990)  and  was  much  further  developed 
by  Sen  and  Stoffa  (1991).  The  Very  Fast  Simulated  Annealing  (VFSA)  modification  of 
Ingber  (1989)  accelerates  the  conversion  by  progressively  focusing  the  model  search  space 
onto  the  better  fitting  solutions.  The  SA  algorithm  used  here  is  an  extension  of  all  of  these 
previous  methods. 

The  SA  algorithm  used  here  has  the  option  of  using  a  different  temperature  (or 
convergence  criteria)  for  each  layer  in  the  earth  model.  By  using  lower  temperatures  for  the 
upper  layers,  the  model  parameters  for  these  layers  can  be  fit  prior  to  fitting  the  lower 
layers.  This  is  similar  to  the  technique  of  layer  stripping.  The  upper  layers,  especially  the 
sea  floor,  usually  give  stronger  reflections  than  lower  layers.  Fitting  the  stronger  features 
from  the  upper  layers  is  usually  easy  and  can  be  done  quickly.  Then  with  the  solution  for 
the  upper  stracture,  the  lower  model  parameters  are  easier  to  fit.  The  second  sigmficant 
new  option  is  for  fitting  the  waveform  envelope  rather  than  the  waveforms.  Fitting  the 
waveform  envelopes  greatly  smoothes  the  residual  function  particularly  for  high  frequency 
and  limited  bandwidfii  signals  making  the  inversion  much  faster  (Wood  and  Lindwall, 
1996). 


Demonstrations 

The  first  inversion  demonstration  is  on  synthetic  data  (Figure  la).  Testing  an  inversion 
algorithm  on  synthetic  data  where  the  exact  solution  is  known  tests  the  inversion  method 
rather  than  data  collection  or  processing.  The  environmental  model  used  for  the  synthetic 
data  in  Figure  1  has  15  layers  with  six  variables  for  each  layer;  these  variables  are: 
thickness,  density,  compressional  velocity  (Vp),  shear  velocity  (Vs),  compressional 
attenuation,  and  shear  attenuation.  The  search  was  done  only  over  Vp  and  layer  thickness 
for  the  upper  3  sediment  layers  in  order  to  reduce  computation  time.  All  other  model 
variables  were  set  to  the  trae  values.  Within  a  Hmited  time  (877  iterations)  the  inversion 
found  a  very  good  solution  shown  in  Figure  Id.  The  residual  amplitude  (Figure  Ic)  was 
less  than  ten  percent  of  the  data  amplitude.  This  first  inversion  fit  the  waveform  rather  than 
the  envelope  as  was  done  for  the  two  field  data  cases. 

The  second  inversion  demonstration  is  with  field  data  (Figure  2a)  from  a  750  m  water 
depth  site  in  the  Eel  River  -  STRATAFORM  region  offshore  Eureka,  California  using  ship 
deployed  sonobuoys  and  a  watergun  sound  source.  The  watergun  produced  a  very  sharp 
impulse  so  source  deconvolution  was  not  used.  The  data  resolution  was  lowered  with  a  50 
Hz  low  frequency  bandpass  filter  (from  a  2  kHz  sample  rate)  to  allow  for  simpler  models 
and  an  easier  inversion.  Source  to  receiver  offsets  were  carefully  determined  from  the  direct 
wave  arrival  times.  I  fit  the  envelope  rather  than  the  waveform  since  the  forward  modeling 
algorithm  does  not  include  important  parts  of  the  waveform  in  the  calculation,  specifically 
the  sea  surface  multiples  and  interbed  multiples  from  the  detailed  structures  near  each  of  the 
major  interfaces.  Fitting  the  waveform  envelope  allows  for  a  much  faster  inversion  with 
nearly  the  same  precision  as  a  full  waveform  fit  (Wood  and  Lindwall,  1996).  The  residual 
from  the  inversion  solution  shown  used  5323  iterations  and  has  55  percent  of  the  amphtude 
of  a  random  solution  which  is  a  good  result  for  field  data. 

There  are  only  two  features  in  the  data  (Figure  2a)  that  were  modeled,  the  sea  floor  and 
a  subsurface  reflection  about  0.2  seconds  after  the  sea  floor  reflection.  The  sea  floor  here  is 
755  m  deep,  the  sea  floor  reflection  hyperbola  fit  a  water  velocity  of  1.498  km/s,  the  Vp  of 
the  upper  sediments  is  1.538  km/s.  The  subsurface  reflection  is  0.147  km  below  the  sea 
floor  and  the  sediments  below  have  a  Vp  of  1.558  km/s.  The  Vp  of  the  lower  sediment  is 
poorly  constrained  since  only  the  relative  amplitude  of  the  reflection  is  known  which  is  also 
dependent  on  the  density  contrast,  another  unknown  value.  Velocity  models  from  nearby 
multichannel  seismic  lines  (C.  Fulthorpe,  personal  communication;  J  Yun,  persond 
communication)  have  similar  subbottom  velocities  (Figure  2d).  The  MCS  data  are  adjusted 
so  that  the  sea  floor  depth  is  the  same  as  my  SA  inversion  solution.  The  sediment  Vp  is 
equal  to  the  lowest  value  obtained  from  the  MCS  data  and  aU  but  one  of  the  MCS  shows  a 
significant  reflector  within  50  m  depth  of  the  reflection  fit  by  the  SA  inversion.  Most  of  the 
MCS  data  used  a  168  channel,  2.5  km  long  streamer  giving  much  higher  quality  data  than 
the  sonobuoy  data  used  for  the  SA  inversion.  The  reliability  of  the  SA  inversion  is 
indicated  by  the  overall  agreement  with  the  MCS  data. 

The  third  set  of  data  for  this  demonstration  was  collected  at  a  250  m  water  depth  site  in 
the  Eel  River-STRATAFORM  region  near  Eureka  California  using  ship  deployed 
sonobuoys  and  an  airgun  source.  The  airgun  has  a  lower  frequency  and  narrower 
bandwidth  than  the  watergun  used  at  the  750  m  water  depth  site.  Source  to  receiver  offsets 
were  initially  determined  from  the  direct  wave  arrival  times  but  were  not  accurate  enough 
for  a  good  inversion.  I  had  to  include  the  offsets  as  inversion  variables  adding  twenty 
variables  to  the  existing  ten  model  variables.  Solving  for  thirty  variables  instead  of  ten 
means  not  only  that  each  sweep  (once  through  all  variables)  takes  three  times  as  many 
calculations  but  the  cooling  process  must  also  be  done  much  more  slowly.  There  have  been 
no  studies  of  how  many  more  iterations  are  needed  for  each  new  variable  added  to  an 
inversion  but  my  experience  suggests  that  it  is  a  power  law  function  rather  than  linear.  This 


third  inversion  case  used  35421  iterations  for  thirty  variables  while  the  second  case  used 
5323  iterations  for  seven  variables.  I  again  fit  the  envelope  rather  than  the  waveform.  The 
lower  frequency  of  the  airgun  signal  reduced  some  of  the  complications  of  having  a  limited 
bandwidth  and  the  lack  of  surface  reflections  in  the  synthetic  calculations.  The  shallower 
depth  here  however  complicated  the  analysis  since  the  sea  floor  and  subbottom  reflections 
were  perilously  close  to  the  direct  arrivals  which  saturated  the  sonobuoy  receivers. 
Sonobuoy  records  do  not  have  the  correct  amplitudes  or  waveforms  when  they  are 
saturated  and  the  sonobuoy  electronics  take  several  tenths  of  a  second  to  recover.  The 
distorted  amplitudes  of  the  sea  floor  reflections  and  the  overlapping,  decaying  direct  signal 
make  reliable  inversions  difficult. 

Features  in  the  data  (Figure  3  a)  modeled  include  a  complex  sea  floor  reflection  and  a 
subsurface  reflection  about  0.04  seconds  after  the  sea  floor  reflection.  The  sea  floor  here  is 
245  m  deep,  the  sea  floor  reflection  hyperbola  fit  a  water  velocity  of  1.482  km/s,  the  Vp  of 
the  first  9  meters  of  sediment  is  1.573  km/s  and  then  increases  to  1.612  km/s  in  a  43  m 
thick  layer.  The  Vp  below  this  was  modeled  at  1.618  km/s  but  is  poorly  constrained  for  the 
same  reasons  as  in  the  750  water  depth  site.  Velocity  models  from  nearby  multichaimel 
seismic  lines  (C.  Fulthorpe,  personal  communication;  J  Yun,  personal  communication) 
have  similar  subbottom  velocities  (Figure  3d).  The  MCS  data  are  adjusted  so  that  the  sea 
floor  depth  is  the  same  as  my  SA  inversion  solution.  My  sediment  Vp  values  are  in  the 
middle  of  the  MCS  values.  All  of  the  MCS  data  shows  a  significant  interface  at  the  400-440 
m  depth  range  but  the  reflection  arrival  times  are  the  same  as  the  sea  floor  reflection 
multiple  in  my  data  so  it  is  not  visible  in  this  sonobuoy  data.  Most  of  the  MCS  data  used  a 
168  channel,  2.5  km  long  streamer  giving  much  higher  quality  data  than  the  sonobuoy  data 
used  for  the  SA  inversion.  The  reliability  of  the  SA  inversion  is  indicated  by  the  overall 
agreement  with  the  MCS  data.  A  curve  from  Table  IV  of  Hamilton  (1980)  is  plotted  on  the 
inversion  and  MCS  models  (Figures  2d  and  3d)  to  contrast  Hamilton’s  gradient  with  the 
discontinuous  seismic  models.  Hamilton  fit  a  smooth  polynomial  curve  through  Vp  values 
from  20  sites.  These  two  individual  sites  may  be  approximately  consistent  wiA  Hamilton’s 
Vp  values  but  the  discrete  reflectors  and  the  discontinuous  velocity  -  depth  functions  are 
likely  to  dominate  the  geoacoustic  responses  of  these  two  sites. 


Restrictions  and  Future  Improvements 

Knowing  the  source  function  is  cmcial  for  a  good  waveform  or  envelope  inversion.  Most 
real  sources  (such  as  airguns  and  explosives)  have  waveforms  complicated  enough  to 
obscure  or  mimic  real  interface  reflections.  The  negative  effects  of  a  complex  source 
function  can  be  eliminated  by  either  deconvolving  the  source  function  or  by  including  the 
source  in  the  forward  calculation.  The  SA  inversion  code  described  here  has  only  a  few 
rudimentary  sources  built  into  it,  and  its  utility  would  be  greatly  enhanced  by  improving  the 
source  functions  or  by  building  a  deconvolution  module.  Deconvolution  is  an  inversion 
problem  itself  and  could  either  remove  the  source  and  receiver  responses  from  the  data  or  to 
accurately  determine  the  source  wavelet  for  inclusion  in  the  forward  calculation. 

The  sea  surface  is  a  nearly  perfect  reflector  so  for  shallow  receivers  or  sources,  the 
surface  multiples  will  be  as  strong  as  the  primary.  Inversion  with  a  forward  code  that  does 
not  include  surface  reflections  will  try  to  solve  for  these  multiples  as  part  of  the  sea  floor 
reflection  response  and  will  put  strong  interfaces  where  none  exist.  The  currently  used 
forward  code  (ray  ID)  does  not  include  surface  multiples.  We  plan  to  include  computation 
of  sea  surface  reflections  in  the  next  version  of  this  code. 

Accurate  inversions  require  offsets  to  within  half  a  wavelength  (1  m  for  some  of  the 
data  shown  here)  and  this  must  be  determined  from  the  data  itself  using  the  direct  wave. 
The  direct  wave  is  not  a  simple  pulse  but  a  complex  waveform  that  changes  due  to  the 
nonlinear  response  of  the  sonobuoys  from  being  saturated  due  to  the  overly  strong  signal. 


The  sonobuoys  response  to  the  direct  arrival  is  due  to  the  overdriving  of  the  electronics 
from  a  strong  signi.  The  recovery  of  the  electronics  depends  on  the  signal  strength.  For 
these  reasons  the  direct  wave  response  changes  enormously  with  offset  and  can  not  be 
done  by  picking  the  peak  amplitude  or  the  first  swing  in  amplitude.  A  program  to  determine 
offsets  by  cross  correlating  the  direct  waves  for  a  window  of  similar  offsets  and  moving 
the  offset  window  over  the  whole  data  set  would  be  much  faster  and  more  reliable  than  the 
manual  offset  picking  and  iteration  that  is  used  now.  This  enhancement  also  should  be 
included  in  further  development  of  the  rapid  inversion  technique. 

The  inversion  technique  presented  here  should  be  enhanced  with  the  addition  of  dip  as 
a  model  variable.  This  would  accommodate  smooth  2-D  earth  models,  models  that  are 
much  more  realistic  for  littoral  regions.  Inclusion  of  dip  is  straightforward  in  ray  trace 
forward  codes  and  the  ability  to  ran  a  line  at  any  azimuth  would  reduce  the  logistics  of  the 
field  activities.  The  data  used  in  this  demonstration  was  acquired  along  strike  (parallel  to 
bathymetry  lines),  a  restriction  that  need  not  remain  in  future  version  of  this  inversion  code. 
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Figure  Captions 


1.  Synthetic  data  calculated  from  a  realistic  marine  sediment  model  (a)  is  inverted  using 
simulated  annealing  (SA)  to  find  a  good  fit  with  limited  computational  resources  (877 
iterations  and  110  minutes  on  a  Sun  SPARC  10).  The  solution  synthetic  is  in  panel  b  for 
comparison.  The  residual  (a  -  b)  is  shown  in  panel  c  and  has  only  9.6%  of  the  amphtude 
of  the  data  (summed  over  the  time  window)  demonstrating  that  the  fit  is  very  good.  Panel  d 
compares  the  known  physical  model  (red  line)  to  the  inversion  solution  (blue  line).  The 
velocity  search  window  was  from  1.5  to  2.5  km/s  and  the  layer  thickness  window  was  ± 
33%  of  the  true  value. 

2.  Field  data  collected  in  the  Eel  River  STRATAFORM  area  at  a  water  depth  of  756  m  (a) 
is  inverted  using  simulated  annealing  (SA)  to  find  a  good  fit  with  hmited  computational 
resources  (5323  iterations  and  27  minutes  on  a  Sun  SPARC  10).  The  solution  synthetic  is 
in  panel  b  for  comparison.  The  residual  (a  -  b)  is  shown  in  panel  c  and  has  55%  of  the 
amphtude  of  a  random  solution.  Panel  d  compares  the  inversion  solution  (thick  line)  to 
velocities  determined  from  several  high  resolution  multichannel  seismic  (MCS)  data  (thin 
fines)  collected  at  nearby  sites  and  adjusted  so  that  the  sea  floor  depths  correspond.  The  SA 
inversion  solution  is  identical  to  the  lowest  MCS  solution  (the  two  fines  are  superimposed 
down  to  905  m)  and  3  MCS  sites  have  a  reflector  between  900  and  920  m  depth.  This 
approximate  agreement  with  the  MCS  data  indicates  that  the  SA  inversion  gives  reliable 
results. 

3.  Field  data  collected  in  the  Eel  River  STRATAFORM  area  at  a  water  depth  of  245  m  (a) 
is  inverted  using  simulated  annealing  (SA)  to  find  a  good  fit  with  limited  computational 
resources  (35421  iterations  and  9.6  hours  on  a  Sun  SPARC  10).  The  solution  synthetic  is 
in  panel  b  for  comparison.  The  residual  (a  -  b)  is  shown  in  panel  c  and  has  80%  of  the 
amplitude  of  a  random  solution.  Panel  d  compares  the  inversion  solution  (thick  line)  to 
velocities  determined  fi'om  several  high  resolution  multichannel  seismic  (MCS)  data  (thin 
fines)  collected  nearby  and  adjusted  so  that  the  sea  floor  depths  correspond.  The  Vp  values 
for  the  upper  150  m  of  sediment  from  the  SA  inversion  are  in  the  middle  of  the  MCS 
derived  values  indicating  a  reliable  result  firom  the  inversion.  The  prominent  reflector  in  the 
MCS  data  at  a  depth  range  of  380-440  m  is  obscured  in  the  sonobuoy  data  by  the  sea 
surface  multiple  so  is  not  modeled  in  the  SA  inversion. 
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